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Abstract 

Recent lattice model calculations have suggested that a full-layered crys- 
tal surface may undergo, under canonical (particle-conserving) conditions, a 
preroughening-driven two-dimensional phase separation into two disordered 
flat (DOF) regions, of opposite order parameter. We have carried out ex- 
tensive classical molecular dynamics (MD) simulations of the Lennard-Jones 
fcc(lll) surface, to check whether these predictions are relevant or not for a 
realistic continuous system. Very long simulation times, a grid of tempera- 
tures from (2/3)T m to T m , and unusually large system sizes are employed to 
ensure full equilibrium and good statistics. By examining layer-by-layer occu- 
pancies, height fluctuations, sublattice order parameter and X-ray structure 
factors, we find a clear anomaly at ~ 0.83 T m . The anomaly is distinct from 
roughening (whose incipiency is also detected at ~ 0.94 T m ), and is seen to 
be consistent with the preroughening plus phase separation scenario. 
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The surfaces of a rare-gas solid such as Ar are reasonably well modeled by those of a 
(truncated) Lennard- Jones (LJ) fee solid. The behavior of the LJ surfaces as a function 
of temperature, particularly in the vicinity of the melting point T m , has been the subject 
of a large number of studies, mostly by classical Molecular Dynamics (MD) Based 
on these studies, the general consensus until recently was that the LJ surfaces begin to 
disorder, with increasing temperature T, in a very gradual manner. As T > (2/3) T m , 
surface anharmonicities B and defects [[!]] first build up. In particular, there is a progressive 
growth in the number of surface adatoms/vacancies in the region T ~ 0.7-^0.8 T m and above. 
Their presence undermines the surface crystalline state ||. Eventually, a microscopic quasi- 
liquid surface film is formed, which leads to surface melting as T m is approached [|J. En 
route to surface melting, a surface roughening transition at Tr < T m should also appear || , 
owing to a step free energy softening in presence of the quasi-liquid film. On Ar(lll), for 
example, it is known that T R ~ 0.94 T m f6j. 

Subsequently, however, experimental evidence has appeared [|7|-[TI|, followed by theoret- 
ical work ||lTl-|T6H which upsets this gradual picture, and suggests that surface disordering 
occurs instead through a singularity. In connection with statistical mechanics models intro- 
duced by Den Nijs in the 80's |I^J, this singularity is probably best understood as "prerough- 
ening" (PR). At PR, the originally ordered flat surface turns into a disordered flat (DOF) 
state, where surface steps proliferate, albeit with strict up-down alternation. In this way 
surface flatness is preserved while still gaining entropy from disorder, in particular from step 
meandering \T%. Due to the resulting checkerboard texture of steps, a DOF surface phase 
exhibits, at least in the lattice models, a striking half-integer occupancy in the topmost layer. 

More recently, grand canonical Monte Carlo simulations [0,0 have confirmed the oc- 
currence of a PR phase transition, and of an associated coverage jump at the LJ(lll) surface, 
in the neighborhood of 0.83 T m . 

The questions we address in this note are the following. Can we first of all detect 
preroughening and roughening in a standard canonical (that is, particle-conserving) realistic 
surface MD simulation? And, if the PR scenario is indeed correct for the LJ surfaces, how do 
PR and the appearance of a DOF phase exactly manifest themselves? Finally, why did such 
an important singularity go unnoticed in so many previous, good-quality MD simulations? 

We surmised recently |20] that the answer might be that in MD, the singularity is masked 
by particle conservation, turning a sharp critical onset into a subtler, Ising-like phase sepa- 
ration. A Monte Carlo study of a lattice RSOS (restricted solid-on-solid) model did indeed 
show that under canonical, particle-conserving conditions, an initially full surface monolayer 
spontaneously phase separates above T PR into two DOF regions, each of which has essen- 
tially population 1/2 in the top layer [2C]. However, lattice models may be oversimplified, 
and a more realistic study of this question is strongly needed. 

Here, we describe new extensive canonical MD simulations of the LJ fcc(lll) surface, 
specifically aimed at understanding whether the DOF phase separation is real or not. An- 
ticipating our conclusions, we shall find that the answer is affirmative. However we also 
find that the evidence can only be obtained by employing large sizes, long simulation times, 
and size scaling, usually not considered in previous work. A corollary is that these kinds of 
precautions, aimed at detecting possible DOF surface phase separation, appear mandatory 
for future studies of warm surfaces using canonical MD simulation. 

We conducted all our simulations in a slab geometry, with N L fully mobile layers of 
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M particles/layer and (x,y) periodic boundary conditions. Three rigid fcc(lll) layers are 
added at the bottom of the slab. To account for thermal expansion, the (x, y) simulation cell 
size was expanded with temperature according to an expansion coefficient extracted from 
an independent set of bulk simulations. Since we are aiming at describing a solid-vapor 
interface at full thermodynamic equilibrium, our simulation allowed an approximately equal 
volume of L J gas above the solid surface. The gas particles were contained and backscattered 
by a perfectly reflecting wall placed at a distance Ny "layers" in the vacuum above the last 
crystal layer. 

We studied in detail three samples: (i) SF (small full) where M = 120 particles/layer, 
N L = 13, N v ~ 10; (it) LF (large full) where M = 504, N L = 25, N v ~ 15; (in) LH 
(large half) where M = 504, Nl = 24.5, Ny — 15. The interparticle interaction was 12-6 
LJ, truncated at 2.5 a. The bulk melting temperature for this system was determined to 
be kBT m ~ 0.7 e, by observing that at this temperature a relatively large number (« 4) 
of melted layers were stable while the rate of growth of the liquid film was maximal (full 
melting of the slab being prevented by the three bottom rigid layers). 

A standard fifth-order predictor-corrector algorithm was used for the time integration of 
Newton's equation, with a time step of 0.01 LJ units, amounting to 2.15 x 10 _14 s in the 
case of Ar {e/ks = 120 K, a = 3.40 A). Temperature control based on velocity rescaling 
was used to heat or cool the system. However, all the equilibration and data collection runs 
were done in constant energy, microcanonical conditions. 

In each case, the standard procedure adopted was to equilibrate the system for 10 5 steps, 
at each temperature T chosen in the interval ~ 0.7 T m -r- T m . During this time (amounting 
to about 2 ns for Ar) some surface atoms evaporated and others recondensed. However, 
their number (less than 10 in all cases) was negligible in comparison with one monolayer, so 
that the solid underneath was basically running under particle conserving conditions. Since 
surface evolution in these conditions is strictly determined by diffusion, long equilibration 
times are dictated by the requirement that the lateral diffusion length of a surface particle 
should be roughly comparable to half the (x, y) cell size. With a diffusion coefficient of the 
order of 0.2 x 10 _4 cm 2 /s for Ar at T m ||21|| , 2 ns are sufficient to ensure equilibration by 
diffusion. Following equilibration, we ran a further 5 x 10 4 4- 10 5 steps for a good statistical 
average of thermodynamical quantities, which we now proceed to discuss. 

(i) Height fluctuations. The surface height fluctuations were obtained by calculating 

«> 2 = (^I>.-'>) 2 ) 

where hi denotes the ^-coordinates of all surface particles i whose identity and total number 
Ns are defined at each given configuration, along with the average height h = (1/Ns) J2i hi- 
Configurations were chosen at regular time intervals (one every 250 steps for SF, every 500 
steps for LF, LH), and (. . .) denotes average over configurations. Surface particles were 
identified as those not totally or even partially shadowed by other particles, when viewing 
the surface from the gas and parallel to the z axis, each particle being represented as a 
sphere of atomic radius r D = 0.9a. 

A flat, defect free vibrating crystal surface should be characterized by 5h 2 ~ d 2 , where d 
is the interlayer spacing. Proliferation of surface adatoms/ vacancies should give rise to an 
increase of 5h 2 , which for increasing size should tend to a finite value, so long as the surface is 
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flat, i.e. non-rough. Surface roughening should be further signalled by a size-dependent 5h 2 
increase, scaling with size as logM. Below roughening and neglecting vibrations, an ideal 
DOF surface, with its half-occupied topmost layer, should be characterized by 5h 2 = (l/4)<i 2 . 
Conversely, a full-layered surface, phase-separated into two DOF domains with a height 
difference of d between them, would instead exhibit a larger 5h 2 ~ (l/2)d 2 . 

Fig. |l| shows the behavior obtained for 5h 2 . There is a clear and sudden change in the full- 
surface behavior near a breakdown temperature ~ 0.83 T m . The height fluctuation 5h 2 grows 
fast for both SF and LF below this temperature; above, it levels off to a value close to that of 
LH. Moreover, comparison of SF and LF indicates a stronger size-dependence of 5h 2 at the 
breakdown temperature. Both features, in accordance with previous discussions based on the 
FCSOS model | |20|| , strongly suggest that preroughening is taking place at the breakdown 



temperature, Tpp ~ 0.83 X^. A strong size-dependence of 6h 2 reappears again at higher 
temperature, compatible with an estimated roughening temperature, Tr ~ 0.94 T m . Both 
values are in excellent agreement with the experimentally established values of Tpa/T m ~ 
69/84, obtained from reentrant layering of Ar(lll) /graphite ||, and of T R /T m ~ 80/84, for 
roughening ||. This is, to our knowledge, the first determination of these temperatures ever 
obtained by direct molecular dynamics simulation. 

Between Tp R and Tr, the canonical SF and LF surfaces should really consist, according 
to the lattice model, of domains made of two kinds of phase-separated DOF regions. The 
height fluctuations do not show strong evidence for that, because 5h 2 remains close to (1/4)gP 
instead of (l/2)d 2 . We will return to this discrepancy later, when discussing evidence for 
surface melting. 

Next we calculated the layer occupancies, shown in fig. 0, as a crucial indicator of phase 



separation p0| . Due to phase separation, we would expect that (neglecting vacancies in 
the second and deeper layers) the concentration of vacancies in the first layer, n v , and of 
adatoms one layer above, n a , should be n v = N v /M = (1 — cr)/2, n a = N a /M = a/2, with a 
the "Ising magnetization" of this problem. In particular we expect a = 1/2 for SF, LF, and 
a = for LH, whence n v = 1/2, n a = for a single DOF phase in LH, but n v = n a = 1/4 
for two phase-separated DOF's in SF, LF. Hence at the onset of DOF phase separation 
these concentrations, ordinarily growing with temperature, should stabilize around 1/4 and 
roughly stop growing until roughening. 

We analyse first the full-layer results, SF and LF. Fig. ^| shows the evolution of the va- 
cancy and adatom concentration with temperature. Beginning with initially small concen- 
trations, the growth rate with temperature has a kink-like feature near 0.83 T m . Comparing 
SF an LF, the kink becomes more pronounced with increasing size, particularly for adatoms. 
Above the kink, there is a visible tendency to stabilize concentrations at a plateau value 
which, for both n a and n v , lies between 0.2 and 0.25, only slightly smaller than the expected 
phase separation value of 1/4. A possible reason why n a and n v might tend to fall slightly 
below 1/4 could possibly be traced to a boundary effect between one DOF domain and the 
other. At the boundary, two parallel steps occur near one another, and parallel step repul- 
sion could be expected to push them apart, thus reducing in the neighborhood both n a 
and n v below 1/4. At much higher temperatures, n v nicely tends to 1/4, while n a grows 
again in the neighborhood of 0.94 T m , where we believe that roughening is taking place. 
We conclude that the overall behavior of adatom/vacancy concentration in the full-layer 
simulations provides a clear evidence for DOF phase separation above the break at 0.83 T m . 
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Next, we consider the half-layer system LH which, by contrast, shows (Fig. 0) n a and 
n v to remain close to and to 0.5, in full equilibrium, across the whole temperature range, 
with no particular features throughout. While the presence of two built-in antiparallel 
steps represents a trivial separation of ordered flat phases at low temperatures, we can 
conclude from the lack of further evolution of adatom/ vacancy concentrations that only 
the full-layered surfaces show a tendency to spontaneously phase separate above TpR, fully 
consistent with expectations. 

It should be possible to identify some of these features by a direct examination of the 
physical structure of the surfaces generated by the simulation. Fig. |^ shows snapshots of 
the full- and half-covered surface at the end of runs below and above 0.83 T m . In spite of 
the fuzziness necessarily present in such snapshots, the following features can be recognized: 
a) the low-T full surface (top left) is relatively crystalline, although with a large number of 
adatoms and vacancies (this is about twice the average value in fig. [|); b) the high-T full 
surface (top right) is now spread on two layers (white and grey), indicating phase separation; 
c) the low-T half-covered surface (bottom left) is also spread on two layers (gray and black), 
again indicating phase separation; d) the high-T half-covered surface (bottom right) is very 
disordered, indicating surface melting. 

These results, providing a vivid picture of what the actual surface instantaneously looks 
like, are in support of the occurrence of PR between the two temperatures. Moreover, they 
indicate melting of the outermost layer above PR, a result predicted theoretically |16| and 
independently obtained by Grand Canonical Monte Carlo simulations |T9|| . 

As a final piece of evidence, we have calculated the sublattice, or DOF, order parameter, 
defined as 

P - / _L V f> i7rh */ d \ 

Ur 7 

The order parameter is finite for an ordered flat surface; it should drop to zero precisely at 
Tpr. Due to a lack of adatom- vacancy symmetry, it should be again finite, but substantially 
smaller, between TpR and Tr. The drop should be either continuous or abrupt, according to 
whether PR was critical or first order. This quantity has the additional merit of correspond- 
ing roughly to the X-ray antiphase scattering amplitude, to be expected for such a surface 



As seen in fig. P has a drop, for both SF and LF, more visible for LF, near 0.83 T m . 
Comparison of the two indicates a weak but nonzero size dependence, suggesting a probable 
weak first order nature for PR in this surface. 

Conversely, the half-layer surface LH exhibits simply a monotonic increase of P, roughly 
up to Tr, where results for all three systems, SF, LF, LH eventually merge. When interpret- 
ing results for LH below TpR and for SF, LF above TpR, we should keep in mind that the 
order parameter average is not particularly meaningful for a phase-separated surface. We 
expect that, as a rough approximation, the true grand-canonical surface would show values 
of P similar to LF for T < TpR, with a jump to those of LH for T > TpR (dashed line). 

Finally, we wish to return briefly to the open problem of anomalously small height 
fluctuations found earlier in the DOF separated region of systems SF and LF. It seems 
possible that an explanation could be found if the top layer did develop a liquid-like character 
above 0.83 T m . In conjunction with relatively small domain sizes, high surface diffusion 
would lead to a globally much more homogeneous, and thus flatter, surface than could be 
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expected on the bare, rigid phase separation picture. The high mobility scenario would apply 
in particular if, as has been repeatedly suggested, the first surface monolayer, or fraction 
of monolayer, is in reality solid below 0.83 T m but melted above p3l , [19|| . We qualitatively 
inspected the diffusivity of our surface atoms and found it indeed to be higher for LH. A 
fully quantitative investigation, reported separately confirmed the existence of a jump 
of the top layer lateral diffusion coefficient by more than a factor 2 between LF and LH at 
T < Tpr. That result does suggest that the first layer is premelting at the same time as it 
is preroughening, and indirectly also provides an explanation for the global homogeneity of 
the DOF phase separated outer surface layers. 

In summary, we have found that a canonical, particle conserving MD simulation for an 
fcc(lll) LJ surface contains clear indications for a) preroughening at Tpp = 0.83 T m ; b) 
roughening at Tp = 0.94 T m ; c) a DOF phase separation between Tpr and Tr. The latter 
result is of conceptual relevance, since it shows that the continuous, off-lattice degrees of 
freedom of the real system do not destroy, with some qualifications, the essence of DOF 
physics, previously established only in much less realistic lattice models. 

Direct experimental observation of DOF phase separation will not be possible on the 
solid rare gas surfaces, where the very high evaporation rates immediately establish grand 
canonical equilibrium. In this sense our study is academic for that system, for which we 



developed separately a grand canonical Monte Carlo approach |]I8],|19|1 . On the other hand, 
the DOF phase separation described here could in real life be realized on metals, under 
conditions where evaporation is irrelevant. 

Independently of experimental accessibility, there is an important methodological mes- 
sage to be drawn from the present work, which is very pertinent to the simulation community: 
surfaces which are likely to undergo a DOF transition should as a rule be simulated grand 
canonically. Canonical molecular dynamics simulation in particular is very dangerous, as 
it may inadvertently describe an unphysical state of the surface in the whole temperature 
interval between preroughening and roughening. 

We are grateful to Santi Prestipino for enlightening discussions. Work at SISSA by 
F. C. was under European Commission sponsorship, contract ERBCHBGCT940636. We 
acknowledge support from INFM PRA LOTUS, and by MURST. C. S. J. is grateful to the 
ICTP for hospitality. 
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FIGURES 



FIG. 1. Height fluctuations 5h 2 as a function of temperature for the small full (SF, black 
squares), large full (LF, black circles) and large half-occupied (LH, white circles) sample. Note 
the gap between LF and LH at low temperature, closing up around 0.83 T m . Other features are 
discussed in the text. 

FIG. 2. Layer occupancies as a function of temperature. Symbols are as in fig. |l[ Note how 
the vacancy and adatom and concentrations of the full-layer samples (SF, LF) tend to stabilize 
in the proximity of 1/4 between about 0.83 T m and 0.94 T m , consistent with the presence of two 
separated DOF phases in this region. 

FIG. 3. Top view of instantaneous snapshots of the three outer layers for systems with initial 
full- (top) and half-layer (bottom) occupation, at 0.76 T m (left) and 0.85 T m (right). Atoms have 
been colored in black, grey or white according to the layer they belong to, in turn determined from 
their z coordinate and the density profile along z of each sample. 

FIG. 4. DOF parity order parameter as a function of temperature. Symbols are as in fig. [l]. 
Note the drop around 0.83 T m for the full-layer samples (SF, LF), indicating tendency toward a 
two-level structure. A hypothetical grand canonical system would approximately follow the dashed 
line, drastically changing the layer population at Tpr in correspondence with a large decrease of 
the order parameter. 
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